%% Carbon Taxes Around the World: Cooperation, Strategic Interactions, and Spillovers
% IMF Economic Review
% Alessandro Moro and Valerio Nispi Landi
% Replication files
% This file computes the final steady state

tau=tau_clean;
tauz=tauz_clean;
mu= 1;                                                  
muz=1; 

r=iota/beta;
rz=iota/betaz;         


options = optimoptions('fsolve','MaxFunEvals',300000,'MaxIter',30000,'TolFun',1e-15,'Display','none');

Y1=[yHss,pHss,pBss,pBzss,yFzss];                       % Initial guess
%Y1=[yHss,1,pBss,pBss,prBss,prBss,prGss,prGss,yHss];                       % Initial guess

Y1=Y1+0.001;

[Y,~,flag] = fsolve(@(YY) find_steady_final(YY,gamma,gammaz,eta,beta,betaz,alpha1,alpha2,alpha1z,alpha2z,chi,...
      GAM,GAMz,n,delta,iota,zeta,csi,kappaA,mu,muz,nu,tau,tauz,A,r,rz,v,g,gz,Az,zetae,zetaez,kappaAz),Y1,options);
yH=Y(1);
pH=Y(2);
pB=Y(3);
pBz=Y(4);
yFz=Y(5);

pF=((1-(1-gamma)*pH^(1-eta))/gamma)^(1/(1-eta));
s=(gammaz*pH^(1-eta)+(1-gammaz)*pF^(1-eta))^(1/(1-eta));
dH=(1-beta/betaz)/GAM;
dF=s/GAMz*(betaz/beta-1);
b=-(v+dH+dF*(1-n)/n);
gdp=pH*yH;
rkB=iota/beta-(1-delta);            % brown return of AE domestic capital
rkB_zd=iota/beta-(1-delta);         % brown return of EME foreign capital
rkB_zu=iota/betaz-(1-delta);       % brown return of AE foreign capital
rkB_zz=iota/betaz-(1-delta);       % brown return of EME domestic capital
rkG=iota/beta-(1-delta);            % green return of AE domestic capital
rkG_zz=iota/betaz-(1-delta);       % green return of EME domestic capital
rkG_zd=iota/beta-(1-delta);         % green return of EME foreign capital
rkG_zu=iota/betaz-(1-delta);       % green return of AE foreign capital
yB=zeta*(pB/pH)^(-csi)*yH;
pB2=pB-tau*(1-mu)*zetae-kappaA/(1+nu)*mu^(1+nu);
pB2z=pBz-tauz*(1-muz)*zetaez-kappaAz/(1+nu)*muz^(1+nu);
pGY=pH*yH-pB*yB;
kB=alpha1*iota/rkB*pB2*yB;
kG=alpha1*iota/rkG*pGY;
i=(kB+kG)*(1-(1-delta)/iota);
kB_zu=alpha2*iota/rkB_zu*pB2*yB;
kG_zu=alpha2*iota/rkG_zu*pGY;
i_zu=(kB_zu+kG_zu)*(1-(1-delta)/iota);
hB=(yB/(A*(kB/iota)^alpha1*(kB_zu/iota)^alpha2))^(1/(1-alpha1-alpha2));
w=(1-alpha1-alpha2)*pB2*yB/hB;
hG=(1-alpha1-alpha2)*pGY/w;
h=hG+hB;
yG=A*(kG/iota)^alpha1*(kG_zu/iota)^alpha2*(hG)^(1-alpha1-alpha2);
pG=pGY/yG;
yBz=zeta*(s*pBz/pF)^(-csi)*yFz;
gdpz=pF*yFz/s;
pGYz=gdpz-pBz*yBz;
kB_zz=alpha1z*iota*pB2z*yBz/rkB_zz;
kG_zz=alpha1z*iota*pGYz/rkG_zz;
kB_zd=alpha2z*iota*pB2z*yBz/rkB_zd;
kG_zd=alpha2z*iota*pGYz/rkG_zd;
i_zd=(kB_zd+kG_zd)*(1-(1-delta)/iota);
i_zz=(kB_zz+kG_zz)*(1-(1-delta)/iota);
c=-(i+pH*g+kappaA/(1+nu)*mu^(1+nu)*yB+b+(1-n)/n*s*i_zd-s*(1-n)/n*(rkB_zd*kB_zd    +rkG_zd*kG_zd)    /iota+(rkB_zu*kB_zu    +rkG_zu*kG_zu)    /iota+rz            *(dH+v)    /iota-r/iota    *(b+    dH    +v)-gdp);
cz=n/((1-n)*gammaz)*(pH/s)^(eta)*(yH-((1-gamma)*pH^(-eta)*(c+i+kappaA/(1+nu)*mu^(1+nu)*yB+i_zu)+g))-(i_zz+i_zd+kappaAz/(1+nu)*muz^(1+nu)*yBz);
wz=(cz^(1/chi)*(1-alpha1z-alpha2z)*(pB2z*yBz+pGYz))^(chi/(1+chi));
hz=(wz/cz)^(1/chi);
hBz=(1-alpha1z-alpha2z)*pB2z*yBz/wz;
hGz=hz-hBz;
pGz=pGYz/yGz;   
e=zetae*(1-mu)*yB;                                                                                                 
ez=zetaez*(1-muz)*yBz; 


yHend=yH;
pHend=pH;
pBend=pB;
pBzend=pBz;
yFzend=yFz;
bend=b;
kBend=kB;
rkGend=rkG;
rkG_zzend=rkG_zz;
tauend=tau;
tauzend=tauz;
iend=i;
i_zzend=i_zz;
cend=c;
czend=cz;



% Save parameters
save paropt yHend pHend pBend pBzend yFzend bend kBend rkGend rkG_zzend tauend tauzend iend i_zzend cend czend

